% Simulation of the two-productivity model in Section 5 of "Repricing
% Avalanches"
% Needs subroutines setparams.m and get_equilibrium.m
% October 2022
% Makoto Nirei

clearvars;

% Set parameter values by setparams.m
setparams;

piv=[0.01 0.03 0.05 0.07 0.09 0.12];   % inflation levels
pilength=length(piv);

% Calvo timings. The same random series is reused for various levels of pi
Fcalvo=rand(Z,1);
calvot = log( 1-Fcalvo )/(-mu*n); % Elapsed time between two Calvo shocks
calvot = cumsum(calvot); % Calendar time of Calvo shocks
maxmonth=floor(calvot(Z)/month_time);   % # of months in the simulation

% Calvo and productivity shocks. Reused for various pi's
ind_c(:,1) = ceil(rand(Z,1)*n); % index of Calvo-hit firms
ind_c(:,2) = (rand(Z,1)<zeta); % productivity shock. 0=non-switch, 1=switch

calvom=zeros(maxmonth,1);   % calvom(i): cumulative # of Calvo shocks by i-th month
calvom(1)=sum(calvot<month_time); % calvom(1)= # Calvo in the 1st month
for i=2:maxmonth
   calvom(i) = sum(calvot<(month_time*i));
end

% main loop for each inflation level piv(i_pi)
dlogP=zeros(Z,pilength); nump=zeros(Z,pilength); NpY=zeros(Z,pilength); % preallocation
negnump=zeros(1,pilength); ss_pi=cell(pilength,1);
for i_pi=1:pilength
    tic;
    pi=piv(i_pi);
    
    get_equilibrium; % Compute steady state
    ss_pi{i_pi}=ss;
    
    disp('eta rho zeta mu delta w');
    disp([eta rho zeta mu delta w]);
    disp('a(i) pbar(i) pstar(i) q(i) pu(i)'); % a, \underline{p}_a, p^*_a, \bar{p}_a
    disp([a(1) pb(1) pt(1) q(1) pu(1); a(2) pb(2) pt(2) q(2) pu(2)]);
    disp(pi);
    
    % stationary distribution g(s), using varsigma=0.5
    Gs = @(s,mu,pi,qL,qH) ((qH.^(s*mu/pi)-1)/(qH^(mu/pi)-1) + (qL.^(s*mu/pi)-1)/(qL^(mu/pi)-1))/2;
    
    % shorthands
    ometa=1-eta; % one-minus-eta
    obometa=1/ometa;
    epb=(pb').^ometa; % "ep" = p^{1-\eta}. "pb"=p underbar. "pt"=p top=p*. 
    ept= (pt').^ometa;
    eq=(q').^ometa;
    epu=(pu').^ometa; % "pu"=\bar{p}
    equ=(qu').^ometa; % "eq"=q^{1-\eta}   
    l_pb=log(pb);
    
    psimu_tmp=zeros(n/2,2); psimu=zeros(n,2); epsimu=zeros(n,2);
    % draw initial price profile for each a(i)
    % Using G(s|a) = (q_a^{(mu/pi)s}-1)/((mu/pi)log(q_a))
    % Letting u=rand=G(s|a), obtain p= pb_a q_a^s
    for i=1:2
        psimu_tmp(:,i)=((rand(n/2,1))*(q(i)^(mu/pi)-1)+1).^(pi/mu)*pb(i);
    end
    psimu_tmp = psimu_tmp / mean(mean(psimu_tmp.^(1-eta)))^(1/(1-eta)); % normalize to P=1
    psimu_init=psimu_tmp(:);
    psimu(:,1)=psimu_init;
    epsimu(:,1)=psimu(:,1).^ometa;  % epsimu(:,1) = p_i^{1-\eta}
    epsimu(1:n/2,2)=1;              % epsimu(:,2) = a_i
    epsimu(n/2+1:end,2)=2;          % Initially, the first half firms are a_L, the second half a_H.
    epsold=epsimu(:,1);    
    for z=1:Z % z-th Calvo shock
        ePold = mean(epsimu(:,1));  % aggregate P^{1-eta}(t-)
        epsold(:,1) = epsimu(:,1);
        cfirm=ind_c(z,1);           % index of z-th Calvo-hit firm
        if ind_c(z,2)==0            % No switch of productivity
            epsimu(cfirm,1) = ept(epsimu(cfirm,2)); % reprice to p^*(a)
            ePnew = mean(epsimu(:,1));  % update eP = P^{1-eta}
            deP = ePnew/ePold;      % shift in eP = P^{1-\eta} from eP(-)
            ind2 = find(epsimu(:,1)/deP > epb(epsimu(:,2))); % update ep and find firms hitting threshold.
                                    % note that eP shifts down when eP shifts up since 1-eta<0 
            numf = 1;   % number of repricing firms including Calvo
            epsimu(ind2,1)=epsimu(ind2,1).*eq(epsimu(ind2,2)); % "ind2" firms reprice to pq
            flag = length(ind2);    % # of ind2 firms
            while flag>0 && numf < n    % Avalanche continues until no ind2 firms
                numf=numf+flag;         % accumulate # of repricing firms
                ePnew = mean(epsimu(:,1));  % Update eP
                deP_tmp=ePnew/ePold;        % shift in eP from eP(t-)
                ind2 = find(epsimu(:,1) > deP_tmp * epb(epsimu(:,2)));  % find the next round ind2
                epsimu(ind2,1) = epsimu(ind2,1).*eq(epsimu(ind2,2));    % update ep
                flag=length(ind2);      % # of ind2 firms
            end
        else    % productivity switch
            if epsimu(cfirm,2)==1    % Shock hits a_L firm. Positive tech shock; price cut
                epsimu(cfirm,2)=2;  % Update to a_H
                epsimu(cfirm,1)=ept(2); % Reprice to p^*(a_H)
                ePnew = mean(epsimu(:,1));
                numf = -1;          % Negative avalanche = price cuts
                deP = ePnew/ePold;
                ind2 = find(epsimu(:,1)/deP < epu(epsimu(:,2)));    % Upper threshold for p; Lower threshold for ep
                epsimu(ind2,1) = epsimu(ind2,1).*equ(epsimu(ind2,2));   % Update p using pu
                flag = length(ind2);                                % number of price cutting firms
                while flag ~= 0 && -numf < n
                    numf=numf - flag;                               % price cut is stored as a negative numf
                    ePnew = mean(epsimu(:,1));
                    deP_tmp=ePnew/ePold;
                    ind2 = find(epsimu(:,1) < deP_tmp * epu(epsimu(:,2)));  % using threshold pu
                    flag=length(ind2);
                    epsimu(ind2,1) = epsimu(ind2,1).*equ(epsimu(ind2,2));
                end
            else                    % Shock hits a_H firm. Negative shock. Price increase.
                epsimu(cfirm,2)=1;  % Update to a_L
                epsimu(cfirm,1)=ept(1); % Reprice to p^*(a_L)
                ePnew = mean(epsimu(:,1));
                numf = 1;
                deP = ePnew/ePold;
                ind2= find(epsimu(:,1)/deP > epb(epsimu(:,2)));
                epsimu(ind2,1)=epsimu(ind2,1).*eq(epsimu(ind2,2));
                flag=length(ind2);
                while flag>0 && numf < n
                    numf=numf+flag;
                    ePnew = mean(epsimu(:,1));
                    deP_tmp= ePnew/ePold;
                    ind2 = find(epsimu(:,1) > deP_tmp * epb(epsimu(:,2)));
                    epsimu(ind2,1) = epsimu(ind2,1).*eq(epsimu(ind2,2));
                    flag=length(ind2);
                end
            end
        end        
        deP = ePnew/ePold;
        epsimu(:,1) = epsimu(:,1)/deP;      % Update the price profile for the next Calvo event
        dlogP(z,i_pi) = log(deP)*obometa;   % Store the price increase at this instant
        nump(z,i_pi)=numf;                  % Store # of repciring firms at this instant including Calvo-hit firm
        NpY(z,i_pi)= mean(epsimu(:,1).^(-eta/(1-eta)) ./a(epsimu(:,2))' );  % N/Y
    end
    negnump(i_pi)=sum(nump(nump(:,i_pi)<0,i_pi))/sum(abs(nump(:,i_pi)));    % Fraction of price cuts
    toc;
end
%save('simu.mat','-v7.3');